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SUMMARY 

The vortex-scalar element method, a scheme which utilizes vortex elements 
to discretize the region of high vortlclty and scalar elements to represent 
species or temperature fields, Is utilized In the numerical simulations of a 
two-dimensional reacting mixing layer. Computations are performed for a diffu- 
sion flame at high Reynolds and Peclet numbers without resorting to turbulence 
models. In the nonreacting flow, the mean and fluctuation profiles of a con- 
served scalar show good agreement with experimental measurements. Results for 
the reacting flow Indicate that for temperature-independent kinetics, the chem- 
ical reaction begins Immediately downstream of the splitter plate where mixing 
starts. Results for the reacting flow with Arrhenius kinetics show an Ignition 
delay, which depends on the reactants temperature, before significant chemical 
reaction occurs. Harmonic forcing changes the structure of the layer, and con 
comltantly the rates of mixing and reaction, In accordance with experimental 
results. Strong stretch within the braids In the nonequilibrium kinetics case 
causes local flame quenching due to the temperature drop associated with the 
large convective fluxes. 


1. INTRODUCTION 

Turbulent diffusion flames have been the subject of extensive experimental 
and theoretical Investigations during recent years (for a review, see Bllger 
(ref. 1)). In most of the theoretical work, turbulence models are used to 

close a system of averaged transport equations which describes the statistical 

behavior of the aerothermodynamlcal variables. Moment methods (ref. 2), eddy 
break-up and mixing controlled models (ref. 3), flame sheet approximation 
(ref. 4), assumed probability density function (PDF) shape methods (ref. 5), 
solutions based on modeled joint PDF of scalar quantities (refs. 6 and 7), and 

based on modeled joint PDF of scalar and velocity (ref. 8) are examples In 
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which turbulence modeling have been used for the closure of equations governing 
the statistical quantities. Much effort has gone Into constructing accurate 
models and In obtaining results that are In agreement with experimental meas- 
urements. However, In complex systems, modeling Is difficult because of our 
lack of knowledge on the detailed dynamics of the flow. Furthermore, since 
most of the Interesting dynamical behavior of the flow Is modeled a priori, 
such features are not exhibited from the results of numerical computations 
based on turbulence models, and thus can not advance our understanding of tur- 
bulent combustion. 

The progress In numerical methods and the availability of supercomputers 
have had a major Impact on turbulence research. Improved accuracy of the 
numerics and Increased storage and computational speed have made It possible 
to solve the appropriate transport equations governing turbulent combustion 
Indirectly without the need for modeling over some limited parameter range. 

Such nearly model-free "simulations," In comparison with calculations utilizing 
turbulence models, have the advantage that the dominant physics of the problem 
Is not modeled a priori, but Is recovered directly from the computed results. 
Their results can be used to understand many Important mechanisms of turbulent 
transport and Its direct Influence on chemical reactions. Furthermore, since 
the Instantaneous behavior of the variables are known at all points and at all 
times, accurate simulations offer a good method of probing the flow when exper- 
imental techniques may fall. There are, however, some limitations on the range 
of turbulent scales that can be resolved accurately by model free simulations. 
Therefore, there Is a need to validate the results of the simulations by a 
direct comparison with experimental measurements. With such validations, ab 
Initio predictions can ultimately be a reality. 

Numerical methods have been used In a variety of forms for the simulation 
of turbulent flows in complex configurations. A recent survey can be found in 
review articles (refs. 9 and 10). In reacting flow, three approaches are used: 
(1) finite difference methods, (2) spectral methods; and, (3) vortex methods. 

In the first approach, the variables are defined on a grid and the transport 
equations are approximated by discretizing the derivatives on the grid nodes. 
Examples of this approach can be found In the work of Corcos and Sherman (ref. 
11) who used a projection method to study the temporal evolution of a periodic 
shear layer, and In Grlnsteln et al . (ref. 12) who used a flux-corrected trans- 
port scheme to simulate the development of coherent structures in a two- 
dimensional spatially evolving shear layer and examined their effect on mixing. 

In spectral methods, the variables are expanded In series of harmonic 
functions that satisfy the differential equations on a number of collocation 
points. Riley et al. (ref. 13) used a pseudo-spectral scheme to study a three- 
dimensional temporally-evolving reacting mixing layer assuming a constant reac- 
tion rate, constant density, and no heat release. McMurtry et al. (ref. 14) 
considered the effects of the chemical heat release on the fluid dynamics of a 
two-dimensional mixing layer for a constant reaction rate. The interplay 
between fluid dynamics and the chemical reaction is Investigated under these 
conditions. Givi et al. (ref. 15) used the same method to compute a two- 
dimensional mixing layer with an Arrhenius chemical reaction and constant dens- 
ity to assess the effects of large coherent structure on the local extinction 
of the flame. Extension to spatllly-growlng layers was Initiated by Givi and 
Jou (ref. 16) using a hybrid pseudo-spectral second order finite difference 
scheme. In all cases, the Reynolds number was kept at small values, 0(100), 
limited by the grid resolution and the number of harmonic modes. 
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In the third approach, vortex methods are used. These schemes are grid 
free, the transport of the variables takes place In a lagranglan form, and the 
solution Is not restricted by the geometry of the confinement. Therefore they 
can provide accurate simulations for high Reynolds number, spatially growing 
flows. Moreover, vortex methods optimize the computational efforts by distri- 
buting computational elements around regions of high vortlclty. The applica- 
tion of the method In thin premixed flame calculations with a finite density 
jump has been reported by Ghonlem et al. (ref. 17) and Sethlan (ref. 18), among 
others. In these calculations, the vortex method was employed to compute the 
flow field, and the dynamic effect of combustion was represented by the propa- 
gation of a thin Interface at the laminar burning velocity acting as a volumet- 
ric source. 

Vortex methods were also used In simulating diffusion flames In connection 
with a finite-difference approach for the treatment of the scalar variables. 
Ashurst and Barr (ref. 19) used the vortex method to compute the hydrodynamic 
field and an Eulerlan flux-corrected transport algorithm to compute the diffu- 
sion and convection of a conserved Shvab-Zeldovlch scalar approximating the 
shape and convolution of the flame in the limit of Infinitely fast chemical 
reaction. Lin and Pratt (ref. 20) used the random vortex method to simulate 
the large-scale motion and a Monte-Carlo method to calculate the time-dependent 
probability density function of the scalar quantities for both gaseous and 
aqueous mixing layers. The PDF transport equation, however, required a closure 
model for the molecular mixing term. 

From this short review. It Is clear that numerical simulations have played 
an Important role In elucidating the physics of turbulent reacting flows, and 
that there Is a continuing need for more direct simulations In order to explain 
better some of the Interesting physical phenomena that have been observed in 
laboratory experiments. 

In this work, we extend the vortex method to study nonpremlxed chemical 
reactions. A vortex-scalar element method Is developed to treat both the 
hydrodynamic and the scalar field In a Lagranglan sense. The fact that a chem- 
ical reaction Is truly a Lagranglan process, l.e.. It occurs when the particles 
(or macroscopic elements) Interact as they flow, motivate the Implementation 
of Lagranglan methods for simulations of high Reynolds number reacting flows. 
The method Is capable of handling a wide variety of initial and boundary condi- 
tions and Is not limited to simple flow boundaries. In this paper, we concen- 
trate on the formulation of the model and the numerical schemes, and present 
some preliminary validation studies and Interpretations of the results. 

In section 2, the geometrical configuration of a spatially evolving mixing 
layer Is presented, and the formulation of the problem and of the scheme are 
described. Results of some sample calculations are given In section 3. Compu- 
tations of a nonreacting mixing layer Is performed first In order to check on 
the accuracy of the method by comparing Its results with experimental measure- 
ments at the same conditions. Preliminary results of a reacting mixing layer 
simulation In which the two reactants are Introduced In different streams are 
presented next. Both constant rate kinetics and temperature dependent kinetics 
are considered. In both cases, the Influence of the coherent structures on the 
finite rate chemistry Is assessed and In the second case, the nonequilibrium 
effects In the reaction rate are examined. In the constant rate kinetics cal- 
culations, the Influence of harmonic forcing at the Inlet of the mixing layer 
Is Investigated. This study was motivated by recent experimental observations 
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of Roberts and Roshko (ref. 22) and numerical computations of Ghonlem and Ng 
(ref. 23). The paper Is concluded In Section IV with a summary of our new 
results and suggestions for future developments. 


2. FORMULATION AND NUMERICAL SCHEME 


A two-dimensional, confined, planar mixing layer Is considered. A sche- 
matic diagram of the flow field Is shown In figure 1. Two Initially unmixed 
reactants, fuel F and oxidant 0, are present at small concentrations In the top 
high speed stream and bottom low speed stream, respectively. We make the fol- 
lowing assumptions: (1) the heat release Is low so that Its effect on the 

dynamics of the flow Is negligible; (2) the Mach number Is small; (3) the free 
stream concentrations of F and 0 are equal and constant; (4) the molecular d 1 f - 
fuslvltles are equal and constant; (5) the viscosity Is the same In both 
streams; and (6) the chemical reaction between F and 0 Is single step, Irre- 
versible, and second order. The density Is, therefore, constant, and the 
transport equations of the hydrodynamic field and the scalar -- temperature or 
species — fields are decoupled. The equations governing this system are: 


F * 0 


7 \|r = -w(X f t) 


do> 1 2 

at + u * - Ri 7 “ 


U ♦ U . VT = ^ V 2 T + Q Da W 


it 1 * u ' ,C J = pTii v % - Da * 


ac p 12 • 

w + u • * c p - ptu v c p + Da w 


(1) 

( 2 ) 

(3) 

(4) 

(5) 

( 6 ) 


where P Indicates products and W = cpco exp(-Ta/T) Is the reaction rate, 
written In terms of the rate of generation of products per unit mass, 
u = (u,v) Is the velocity, x = (x,y) and x,y are the streamwlse and cross 
stream directions, respectively, t Is time, Is the stream function defined 

such that u = a^/ay and v = -a»|r/ax, w = vxu Is the vortlclty, c Is the 
concentration per unit mass, T Is temperature. V = (a/ax, a/ay), and 
v 2 = a 2 /ax 2 + a 2 /ay 2 . Variables are nondlmenslonallzed with respect to the 
appropriate combination of the total shear AU = U1 - U2, the channel height 
H, the free stream concentration of F, cp 0 , the free stream temperature at 
x = 0, To. In equation (5), j = F or 0 for fuel and oxidizer, respectively. 
Re - AU H/v Is the Reynolds number, where » Is the kinematic viscosity. The 
reaction rate constant k = A exp(-Ta/T) where A Is the frequency factor, and 
Ta Is the activation energy, nondlmenslonallzed with respect to (RTo), R 
being the gas constant. Q Is the enthalpy of reaction, nondlmenslonallzed 
with respect to CpTo, where Cp is the specific heat at constant pressure. 

Pe = AU H/a Is the Peclet number, where a Is the thermal dlffuslvlty. 
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Da = A cp 0 H/AU Is the first Damkohler number. D Is dlffuslvlty and 
Le = a/D Is the Lewis number. 

Since equations (4) to (6) are similar, there Is no need to solve them all 
If the scalar concentrations cp, eg, and cp are normalized In such a way 
that their Initial and boundary conditions are Identical. This Is accomplished 
by the use of Shvab-Zeldovlch transformation (ref. 1). Introducing conserved 
scalars Bpp = cp + cp, and Bgp = 1 - (eg - cp), we get: 



VB 


_J 

j " Pe Le 



(7) 


for j = FP or OP. Since Bpp and Bop have the same Initial and boundary 
conditions, Bpp = Bop « B. The finite rate kinetics effects can be taken Into 
account by considering the transport equation for the product of chemical reac- 
tion, equation (6), and equation (7) for a conserved scalar. If the Lewis num- 
ber Is unity, another conserved scalar can be Introduced, Bpj = cp - T/Q, and 
the solution of equations (6) and (7) for c p and B will determine the behav- 
ior of all the scalar equantltles, cp, eg, cp, and T. 


2.1 The Vortex Scheme 

In the vortex method, the vortlclty field Is represented by a finite num- 
ber of vortex elements of finite cores: 

w (x,t) = i rV« 2 f (x - xi) (8) 

where = J u dA, Is the circulation of a vortex element and A Is the core 
radius, while xi Is the center of the element, f represents the vorltlclty 
distribution associated with a vortex element, or the core function (Chorln 
(ref. 23), Hald (ref. 24), and Beale and Majda (ref. 25).) The velocity field 
Is obtained by solving equation (2) using the discrete vortlclty distribution. 

u = l K(x - xiMx - xi) + Up (9) 

where K(x) « -(y,-x)/r 2 Is the kernal of the Poisson equation, 

*( x) = J r f(r) dr Is the circulation within r, and r = |x|. u p Is an 
Irrotatlonal velocity field added to satisfy the potential boundary condition; 
Up = v<j> where v 2 <t> = 0 and u.n =0 on solid boundaries while u.n = U at 
tne Inlet, n Is the normal unit vector. For the confined shear layer, the 
boundary condition at x = 0 Is: u = U1 for y > 0 and u = U2 at y < 0, 

while y = 0 Is a vortex sheet of strength AU = U1 - U2. 

In this work, we use Ranklne vortex elements, l.e., the vortlclty of an 
element Is constant within the core and zero outside, f(r) = 1/ir for r < 6 
and f(r) = 0 for r > 6. Correspondingly, K(r) = r 2 /2ir for r < i and 
k = 1 for r > S. Moreover, the potential velocity field Is obtained by con- 
formal transformation. Thus, the physical plane Is mapped onto the upper half 

plane and Image vortices are used to satisfy the potential boundary conditions. 
The form of the mapping function for the confined shear layer Is given by 
Ghonlem and Ng (ref. 22). 
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The motion of the vortex elements must be constructed such that the vor- 
tlclty field satisfies equation (3). This Is accomplished by solving this 
equation In two fractional steps: 

do> 

Convection: + u • v« * 0 (10) 

Diffusion: ff = (11) 

In the first step, the convective transport of vortlclty Is Implemented In 
terms of the Lagranglan displacement of the vortex elements using the current 
velocity field computed from equation (9). In the second step, the solution 
of the diffusion equation Is simulated stochastically by the random walk dis- 
placement of the vortex elements according to the appropriate population. 

Thus : 


x^(t + At) = x^(t) + I k u(x 1k )At = (12) 

for 1 = 1,2 N, where Is a k-th order time-integration scheme and 

m Is a tw o-dimen sional Gaussian random variable with zero mean and standard 
variation V 2At/Re. For more details, see Ghonlem and Ng (ref. 22), Ghonlem 
and Gagnon T ( ref. 26). 

The no-slip boundary condition at the walls Is satisfied by generating new 
vortex elements to cancel the Induced velocity by the vortlclty field. Here, 
we generate vortlclty only at the point of separation, l.e., at the tip of the 
splitter plate since the growth of the boundary layers along the channel walls 
at these high Reynolds numbers Is small. At each time step, the new vortlclty 
Ar = -All Urn At, where Urn = (U1 * U2)/2, Is consigned to No elements of 
strength Ar/No and added to the field at points Ax = Um/No apart down- 
stream of x - 0. 

The effect of the numerical parameters on the accuracy of the results was 
Investigated by Ghonlem and Ng (ref. 22). Their results emphasized the Impor- 
tance of using a high order time-integration scheme with k = 2 to avoid 
excessive numerical diffusion In the vortlclty field. The value of No = 6 was 
also found to be appropriate In order to obtain well-defined eddy structures 
after the rollup and the first two pairings. The second pairing is accom 
pllshed within the domain of 0 < x < 6, therefore the computational domain was 
limited to X max = 6. Downstream of X max , the vortlclty was deleted. 

Varying X max showed that the effect of deleting the vortex elements propa- 
gates about one channel highest upstream, hence the results are accurate only 
for 0 < x < 5. 


2.2 THE SCALAR FLEHENT METHOD 

In this scheme, which Is a two-dimensional extension of the random element 
method of Ghonlem and Oppenhelm (ref. 27), the scalar field Is represented by a 
set of elements each carrying a finite amount of the scalar field. 

s(x,t) = z Si i(x - xi ) (13) 
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where s Is a scalar field, being the temperature of species concentration, s^ 
Is the strength of an element, defined as the amount of scalar carried by this 
element and «(.) Is the Dirac delta function, s^ = 1/6A J s(x,t) dA, where 
«A = ixay, and ix and Sy are the distances between the centers of neighbor- 
ing elements In the streamwlse and cross stream directions., respectively, and 
x\ Is the center of the element. If s Is the active scalar. Its transport 
Is governed by: 


ff ♦ u • n - L , 2 S ♦ i (14) 

where Se Is the ratio between the diffusive and convective time scales of 
transport of s, Se = Pe for s = T, and Se = Pe Le If s = c. In the 
scalar element method, this equation Is solved In three fractional steps: 


3 S 

Convection: ^ + u • vs = 0 

(15) 

Diffusion: ff - ^ Vs 

(16) 

Reaction: || = W 

(17) 


Convective and diffusive transport are taken Into account In a similar way 
as In the vortex method, l.e., by the Lagranglan motion of the scalar elements 
using the velocity field u, computed using equation (9), and the random walk 
displacement of the elements using a set of Gaussian random variables with zero 
mean and standard deviation y2At/Se (Ghonlem and Sherman (ref. 28)). If xi 
Is the center of the element 1, then, 

Xi(t + at) = xi(t) + z u(xik) at + ^1 (18) 

Chemical reaction changes the amount of reactants carried by the element 
according to the Integration of equation (17). 

si ( t + at) = si(t) + w at (19) 

However, the reaction occurs only when the element are close enough for molecu- 
lar mixing to affect their composition. Therefore, at every time step, the 

distance between the centers of each two elements of F and 0 

axil * Ixi - xi I Is computed. If axij & «q, where 3 O(iySe) i s the 
diffusion length scale, the composition of each of the two element changes 
according to equation (19). The Initial distance between neighboring elements 
must be small enough to allow enough Interactions between the elements. This 
limits the maximum value of the Peclet number that can be economically used In 
the computations to 0(1000). 

The scheme, while providing an approximate solution of equation (12) In a 
stochastic sense, mimics closely the actual physics of the reaction process. 
This Is achieved by using the Lagranglan formulation of the transport equations 
and dealing with the chemical production terms In Individual particles. 
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3. RESULTS AND DISCUSSION 


The computer code, developed by Ghonlem and Ng (ref. 22) for vortex simu- 
lation of a nonreacting shear layer, was vectorized In order to take advantage 
of the computational capability of a CRAY-XMP. The scheme, being explicit In 
time and requiring mostly nonrecursive computations, can utilize this capabil- 
ity efficiently. The dynamics for the nonreacting layer was Investigated In 
detail In the work of Ghonlem and Ng (ref. 22). Here we concentrate on results 
pertaining to mixing and to chemically-reacting layer. 


3.1 Nonreacting Mixing Layer 

Results of a typical simulation, presented In terms of the velocity and 
location of all vortex elements used In the computations, are shown In 
figures 2 to 4 for the cases of Re = 24 000, Re = 4000, and Re = 1000, respec- 
tively. Each vortex element Is depicted by a point, while Its velocity rela- 
tive to the mean velocity Is represented by a line vector starting at the 
center of the vortex element. The velocity ratio across the layer at the Inlet 
Is U2/U1 = 1/3. 

Results show the formation of large vortex eddies by the rollup of the 
vortlclty layer that emanates at the splitter plate, and the subsequent pair- 
ings of these eddies Into larger structures. The rollup of the shear layer was 
Investigated In Ghonlem and Ng (ref. 22) by analyzing results at a wide range 
of the Reynolds number and at different boundary conditions. Their analysis 
show that: (1) the rollup Is due to the growth of perturbations by the Kelvln- 

Helmholtz Instability mechanism, and the shedding frequency corresponds to the 
most unstable frequency predicted from the linear stability analysis of a spa- 
tially growing layer; (2) pairing, which Is associated with the local subhar- 
monic perturbations, results In a step-wise Increase In the size of the 
vortlclty layer as two eddies merge; (3) The two sources of the subharmonic 
perturbations are the downward motion of the layer and the monotonic growth In 
the size of the eddies downstream; (4) the Intrinsic dynamics of the Instabil- 
ity Is not strongly affected by the value of the Reynolds number, except that 
at the low Reynolds number the eddies are slightly larger due to the dlsperlson 
of vortlclty by diffusion; and (5) the computed velocity statistics show good 
agreements with experimental data, Indicating that the fundamental mechanisms 
of the shear layer are two-dimensional and, hence, the numerical scheme Is cap- 
able of predicting the large scale features accurately. 

The study entrainment, a passive conserved scalar with a normalized con- 
centration value equal to zero In the high speed stream and equal to one In the 
low speed side Is Introduced at the Inlet section. At each time step, 19 ele- 
ments are Introduced In each stream. The Initial distance between two neigh- 
boring elements In the cross stream direction Is taken as 6y = 0.021. The 
time step At = 0.1, thus the distance between the elements In the streamwlse 
direction Is Sx = 0.05 on the average. Since diffusion Is more critical In 
the cross stream direction, 6y Is chosen to be smaller than <5x. A case with 
4y = 0.016, using 25 elements In each stream was computed, showing no signifi- 
cant change In the overall behavior. 

Figures 5 to 7 are obtained for Reynolds number, Peclet number, and veloc- 
ity ratio 10 000, 4 000, and 1/2, respectively. Figures 5 and 6 show the 
velocity and location of all the vortex and scalar elements respectively, while 
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figure 7 exhibits the strength of each of the scalar elements at the nondlmen- 
slonal times of t = 28, 29 and 30. In figure 6, the dots represent the fluid 
from the high speed side with normalized concentration c = 0, and the open 
circles represent the fluid from the low speed side with c = 1. This figure 
Indicates that the rollup of the vortices and their subsequent pairing entrains 
fluid from both sides of the free streams Into the cores of the vortlclty 
layer, which results In the enhancement of mixing between the two streams. 
Entrainment asymmetry Is observed as more fluid from the high speed side Is 
present In the low speed side than the opposite (Koochesfahanl (ref. 29)). 

The Instantaneous profiles of the concentration field are averaged over a 
long-time period and the statistical values are compared with experimental data 
In figures 8 and 9. Figure 8 shows the mean value of the concentration, c m , as 
a function of (y - y 0 )/(x - x 0 ), where y 0 Is measured at c m = 0.5 and x 0 
Is the virtual origin of the mixing layer based on the mean concentration pro- 
file (In the calculation, x 0 =0). In this figure, the solid line Is the 
computed mean concentration at x = 4 and the data points are obtained from 
recent experimental measurement by Masutanl and Bowman (ref. 30) for a dilute 
nonreacting mixing layer with the same velocity ratio. Figure 9 shows a com- 
parison betw een the c omputed and measured mean fluctuations of the concentra- 
tion, c"'2 = (c - 0 ^) 2 . it Is evident from the two figures that both the mean 
and the second moment of the conserved scalar across the width of the shear 
layer are accurately predicted by our computations. 

We note that the results In figures 8 and 9 are In better agreement with 
experimental data than those previously predicted by Givi et al. (ref. 31). 

In these calculations, a k-e turbulence model and a gradient diffusion model 
for turbulent transport of the scalar mean, moment, and probability density 
function was utilized. In the k-e calculations, the concentration fluctua- 
tions exhibit a fairly smooth bell-shaped profile with a much less clear double 
"hump" In the middle region. Indicating poor agreement near the high speed 
stress. The present calculations show the two local maxima In the fluctuation 
profiles that correspond to the location where the gradient of the mean value 
Is highest. The same behavior Is observed by the experimental results of 
Masutanl and Bowman (ref. 30) and Batt (ref. 32). It Is clear that, In accord- 
ance with the findings of Broadwell and Brledenthal (ref. 33), the Intermlt- 
tencv caused by the large coherent structures contributes greatly to the sta- 
tistics of this turbulent flow. 


3.2 Reacting Mixing Layer 

In the calculation of a reacting mixing layer, two reactants F and 0 
are Introduced on both sides of the splitter plate. At x = 0, for y > o, 
cp = 1, and eg = 0, and for y > 0, cq = 1, and cp = 0, while c P = 0. As 
reactants are entrained Into the mixing cores of the layer, they diffuse across 
the original Interface and chemical reaction proceeds. The rollup and pairing 
Increases the original length of the Interface by many folds and allow the 
entrained fluid to diffuse along a larger boundary (Ghonlem et al. (ref. 34)). 
During this process, if the Lagranglan elements utilized to represent the 
Interaction between chemically reacting species are brought close enough so 
that the distance between two neighboring elements Is smaller than the charac- 
teristic diffusion length, they react at the rate defined by equation (17). 
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In figures 10 to 12, we present the velocity, location, and the strength 
of the elements In terms of product concentration for the reacting mixing layer 
with constant rate chemical kinetics and temperature-dependent reaction rate, 
respectively. The amount of the products formed due to chemical reaction Is 
presented by the diameter of the circles In the figures, l.e., larger circles 
Indicate more products. In both cases, Re = 10 000, Pe = 4 000, and 
U2/U1 = 1/3 while Le = 1 . In the constant rate kinetics case, the value of 
the Damkohler number Da = 1 and In the temperature-dependent kinetics Da = 
200, Ta = 10, and Q = 5. Note that In both cases the value of the nondlmen- 
slonal kinetic parameters are low enough so that the effects of heat release 
on the fluid dynamics can be negligible. The stiffness of equation (19) for 
large values of the Damkohler number Imposes a restriction on the time step of 
Integration. In these calculations, we found that at = 0.1 Is sufficiently 
small to accurately Integrate the slow chemistry. 

A comparison between the two figures reveal that under Isothermal condi- 
tions, the products are formed as mixing occurs just downstream of the splitter 
plate, while In the temperature-dependent kinetics calculations, there Is an 
Ignition delay before the reactant reach a temperature high enough to allow any 
significant chemical reaction to occur. Once the reaction begins, the mechan- 
ism of product formation and chemical reaction In both cases are asymptotically 
the same. Increasing the Damkohler number to Da = 400 results In a shorter 
Ignition delay, and preheating the reactants by Increasing the temperature at 
the Inlet to T1 = Q/2 while Da = 200, eliminates the Ignition delay as Indi- 
cated In figures 13 and 14, respectively. 

In order to examine the effects of chemical reaction on the transport of 
species, the concentration statistics In the temperature-independent reaction 
case are presented In figures 15 and 16. These figures correspond to the 
ensemble mean and fluctuations In the bottom-stress species concentration In a 
reacting mixing layer with Da = 1 , U2/U1 = 1/2, Re = 10 000, and Pe = 4000. 

A comparison between figures 15 and 8, and between figures 16 and 9 Indicates 
that near the free stream, the chemistry does not affect the statistical behav- 
ior of the species. Near the reaction zone, however, the mean and the rms 
values of the concentration are lower under reacting conditions, while the 
second hump near the high speed stream side of the rms profile In the nonreact- 
ing layer Is eliminated In the reacting flow due to the local consumption of 
the species by chemical reaction. The same behavior was also observed in the 
experiments of Masutanl and Bowman (ref. 30) In a reacting mixing layer under 
Isothermal conditions. Their results, however, can not be compared quantita- 
tively with the present calculations since the values of the chemical param- 
eters employed In the numerical simulation are substantially lower than those 
of the experiment. 


3.3 Effect of Harmonic Forcing 

The dynamic effect of oscillating the upstream side of the layer was 
studied experimentally by several authors, e.g., Oster and Wygnanskl (ref. 35) 
and Roberts and Roshko (ref. 21) and numerically by Ghonlem and Ng (ref. 22). 
Their results Indicate that In the forced case, eddy Interactions follow four 
stages. In the first stage, the layer rolls up at the harmonic of the forcing 
frequency closest to the most amplified mode. In the second stage, a process 
of accelerated pairings yields a large eddy which Is In tune with the forcing 
frequency. This large resonant eddy appears earlier than It would appear In 
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the case of an unforced layer. In the third stage, paring among resonant 
eddies, which represents a neutrally stable mode. Is disabled and the growth 
of the vortlclty layer Is Impaired for several eddies downstream. In the 
fourth stage, the effect of forcing diminishes and pseudo-random pairing Is 
resumed. Moreover, velocity statistics are affected by forcing, and the sign 
of momentum transfer across the layer Is reversed following pairing. Entrain- 
ment of passive particles was found to be commensurate with the development of 
the vortlclty layer. 

In the recent experiment by Roberts and Roshko (ref. 21), It has been 
observed that periodic forcing has a direct Influence on the outcome of chemi- 
cal reaction across a turbulent shear layer. The results of this experiment 
Indicate that when harmonic forcing Is applied, the mixing rate: (1) Is 

Increased In the Initial stages where the resonant eddy Is forming; (2) Is 
decreased In the Intermediate stage which corresponds to the resonant or 
"frequency-locked" region; and, (3) Is the same as that of the unforced layer 
further downstream. In order to characterize these three regions, the 
Wygnanskl-Oster parameter X w = AU nx/Um^ is utilized, where ft Is the forcing 
frequency (ref. 35). Roberts and Roshko (ref. 21) and Browand and Ho (ref. 36) 
show that the three different regions can be classified according to the local 
value of X w parameter. In region I, X w < 1 , the growth rate Is enhanced. 

In region II, X w > 1 , the frequency-locked region, the growth rate Is 
Inhibited. In region III, the growth rate relaxes to that of the unforced 
layer. 

In order to Investigated this phnomenon computationally, the response of 
the reacting shear layer to the application of low frequency, low amplitude 
perturbations on the upstream side of the shear layer Is computed. Streamwlse 
oscillations are applied on both sides of the layer, hence a pressure preturba- 
tlon Is Imposed without changing the vortlclty field. The streamwlse veloci- 
ties are taken as U1 = 1 + a sin (2irftt), and U2 = a 1)2, where a Is the 
amplitude of forcing. 

The normalized distribution of the product thickness along the mixing 
layer for three cases, ft = 0, 0.5, and 1, Is shown In figure 17. In these 
calculations, a * 0.1, and Re = 4000. The figure Indicates that for ft = 1, 
mixing Is enhanced In the Initial part of the layer, 1 £ x £ 2. The resonant, 
frequency-locked region begins at x = 2 and ends at value x ~ 3. In this 
region, mixing Is reduced and Is less than that of unforced mixing layer. 
Downstream of this region, x a 3, mixing rate resumes Its natural growth and 
reaches asymptotically that of the unforced layer. For lower forcing fre- 
quency, ft = 0.5, the same overall behavior Is observed. In this case, how- 
ever, the results of numerical calculations Indicate that the resonant 
frequency-locked region Is approximately In the range 3 £ 4 x £ 4. A compari- 
son between the range of the frequency locked region calculated here with that 
estimated by Browand and Ho (ref. 36) Is shown on table I. Considering the 
fact that our simulations Ignore the effect of small scale three-dimensional 
turbulence motion, and considering the nonuni versallty of the Browand and Ho's 
curve due to Its Independence to experimental conditions and other Important 
nondlmenslonallzed parameters, this agreement Is encouraging. 
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3.4 Effects of Strain Rate 


It has been shown experimentally by Tsujl (ref. 37), numerically by Llew 
et al. (ref. 38), and analytically by Peters (ref. 39), that the strain rate 
has a major Influence on the flame structure, particularly In nonpremlxed sys- 
tems. In the counter-flow diffusion flame experiments of Tsujl (ref. 37), It 
was observed that Increasing the magnitude of stretch near the flame surface 
results In an Increase of the flow of reactants Into the reaction zone. As a 
result, the chemical reaction Is not able to keep pace with the supply of 
reactants, and the reaction rate Is reduced until local flame quenching occurs. 
The analysis of Peters (ref. 39), which Is based on the method of matched 
asymptotic expansion at large activation energy, shows that the mechanism of 
flame extinction can be addressed by examining the local value of the rate of 
scalar dissipation. This parameter Is viewed by Peters (ref. 39) as the 
Inverse of the diffusion time scale. If the local value of dissipation Is 
Increased beyond a critical limit, the heat conducted away from the diffusion 
flame can not be balanced by the heat produced by the chemical reaction. As a 
result, the maximum value of the temperature decreases, and the reaction even- 
tually ceases. 

By Increasing the number of scalar elements to 38 In each stream while 
decreasing the computational domain to X max = 4, and by preheating the Incom- 
ing reactants to T1 = Q/2 to start the chemical reaction Immediately down- 
stream the splitter plate, we were able to observe this phenomenon. Figures 18 
and 19 show the instantaneous velocity and temperature rise, T - T1, of the 
scalar element at times of t = 1 9 and t = 21 , respectively. In this case, 
the Damkohler number, the normalized enthalpy of reaction, the activation 
energy, and the velocity ratio at the Inlet are 50, 8, 20, and 1/3, respec- 
tively. The cross-stream direction Is enlarged by a factor of 2 for the pur- 
pose of clarity. 

The figures show that the number of scalar elements near the braid, which 
Is the thin link between two neighboring cores. Is only a small portion of the 
total number of elements within the computational domain, which reached more 
than 5100. This Indicates an Instantaneous quenching at the stagnation points 
of the layer. Moreover, the temperature and product concentration In the reac- 
tion zone reach a maximum at the core of the eddies where the vortlclty concen- 
tration Is high, while they reach a minimum at the stagnation point within the 
braid between the neighboring cores where the strain and the scalar gradients 
reach their maximum values. This Is consistent with the results of the pseudo- 
spectral calculations of Givi et al. (ref. 15), and with the experimental 
observations of Tsujl (ref. 37) who showed that the local extinction of diffu- 
sion flame occurs mainly at the regions of high dissipation rate. At these 
regions, the temperature tends to decrease, and If It goes below a critical 
characteristic value, the flame locally extinguishes. 

Quantitative analysis of the effects of stretch on the chemical reaction 
Is rather difficult In the context of present algorithm. This Is due to the 
fact that there are very few scalar elements near the regions of high strain, 
and as shown by Ghonlem et al. (ref. 34), most on the elements tend to be con- 
centrated near the regions with low dissipation. Implementation of a numerical 
scheme based on the transport of the scalar gradients, as In Ghonlem et al. 
(ref. 34) can Improve the accuracy of the analysis substantially, particularly 
those associated with the effects of stretch. In this method, the elements are 
concentrated near the regions of large gradients, or high dissipation, and 
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hence a smaller total number of elements have to be considered. The Implemen- 
tation of this method for the numerical simulation of unpremlxed reacting flows 
Is presently underway to study the effect of strain rate more accurately. 


4. CONCLUSION 

In this work, a numerical scheme based on the transport of computational 
element carrying vortlclty and scalar quantities has been developed to simulate 
a reacting planar, two-stream mixing layer with unmixed reactants. The scheme 
solves the transport equations at high Reynolds and Peclet numbers without 
using models for turbulence closure. A Lagranglan stochastic model Is used to 
Implement the chemical reactions for both constant rate kinetics and variable 
temperature Arrhenius reactions. 

In the nonreacting flow simulations, the calculated statistics of the mix- 
ing of a conserved scalar are In good agreement with experimental data. In 
particular, the numerical results show the presence of two maxima In the fluc- 
tuation profile. In the constant rate reacting flow simulation, the effect of 
chemistry Is to smooth out this curve and produce a single maximum, which 
agrees with the experimental observations. Harmonic forcing enhances the mix- 
ing within the accelerated growth zone of the vortlclty layer, while It Impairs 
the entrainment of the unmixed fluid Into the cores In the resonating region. 

As a result, the numerical simulation Indicates a decrease In the rate of pro- 
duce formation In the frequency-looked region, similar to previous experimental 
findings. 

In the Arrhenius, temperture-dependent kinetics, the mechanism of Ignition 
delay, and the effects of reactants preheating on the decease of the duration 
of this delay Is observed. Also, the nonequilibrium coupling between the sca- 
lar dissipation rate and the flame structure Is revealed as quenching fre- 
quently appears within the braids. To describe this phenomenon more 
accurately, work Is underway to construct a higher order scheme which can pro- 
vide better resolution at the regions of strong strain rates. 
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TABLE i. - 


Frequency locked region 

Q 

Calculated 

Measured (ref. 36) 

0.5 

1.0 

3 4 x < 4 

2 < x -< 3 

2.66 < x < 5.33 
1.33 < x < 2.66 


T 6 
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112 « 0.3335, Re = 24000; TIME = 29.00 





U2 = 0.3333, Re = 24000; TIME = 30.00 


FIGURE 2. - VORTICITY FIELD AT Re = 24000, U2/U1 * 1/3. 



U2 = 0.3333, RE - 10000; TIME = 30.00 
FIGURE 3. - VORTICITY FIELD AT Re = 10000, U2/U1 = 1/3. 




U2 = 0.3333/ Re = 4000; TIME = 30.00 
FIGURE 4. - VORTICITY FIELD AT Re - 4000, U2/U1 = 1/3. 





U2 = 0.5000/ RE = 10000; TIME = 28.00 


•<$§> . i ^f& 






U2 = 0.5000, RE = 10000; TIME = 29.00 
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U2 = 0.5000, RE = 10000; TIME = 30.00 
FIGURE 5. - VORTICITY FIELD AT Re = 10000, U2/U1 = 0.5. 
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U2 = 0.5000, Re * 10000; TIME = 28.00 



U2 - 0.5000, Re - 10000; TIME = 30.00 
FIGURE 7. - CONCENTRATION FIELD AT Re = 10000, PE = 4000, U2/U1 = 0.5. 
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FIGURE 8. - NORMALIZED MEAN CONCENTRATION PROFILE AS A FUNCTION OF THE CROSS-STREAM COORDINATE. 
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U2 = 0.3333, Re = 10000; TIME = 20.00 

FIGURE 11. - PRODUCT CONCENTRATION FIELD, Re = 10000, Pe = AOOO, U2/U1 = 1/3, 
ISOTHERMAL REACTING LAYER. 
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U2 = 0.3333, RE = 10000; TIME = 19.00 



U2 = 0.3333, Re = 10000; TIME = 20.00 

FIGURE 12. - PRODUCT CONCENTRATION FIELD, VARIABLE TEMPERATURE REACTING LAVER. 



U2 = 0.3333, Re = 10000; TIME * 20.00 

FIGURE 13. - PRODUCT CONCENTRATION FIELD, VARIABLE TEMPERATURE REACTING LAYER. 



U2 = 0.3333, RE = 10000; TIME = 20.00 



U2 = 0.3333, Re = 10000; TIME = 21.00 

FIGURE 14. - PRODUCT CONCENTRATION FIELD, VARIABLE TEMPERATURE REACTING LAYER. 
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FIGURE 15. - NORMALIZED MEAN CONCENTRATION 
PROFILE AS A FUNCTION OF THE CROSS- 
STREAM COORDINATE. 
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FIGURE 16. - NORMALIZED RMS CONCENTRATION 
PROFILE AS A FUNCTION OF THE CROSS- 
STREAM COORDINATE. 
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FIGURE 17. - VARIATION OF THE PRODUCT 
THICKNESS VERSUS THE DOWNSTREAM 
DISTANCE. 
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U2 = 0.3333/ RE = 10000; TIME = IS. 00 



U2 = 0.3333/ Re = 10000; TIME = 19.00 
FIGURE 18. - TEMPERATURE FIELD FOR REACTING MIXING LAYER. 



U2 - 0.3333/ RE = 10000; TIME = 21.00 



U2 = 0.3333/ Re = 10000; TIME - 21.00 
FIGURE 19. - TEMPERATURE FIELD FOR REACTING MIXING LAYER. 
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